*
* $Id$
*
* $Log: gsrotm.F,v $
* Revision 1.2  2002/12/02 16:37:45  brun
* Changes from Federico Carminati and Peter Hristov who ported the system
* on the Ithanium processors.It is tested on HP, Sun, and Alpha, everything
* seems to work. The optimisation is switched off in case of gcc2.xx.yyy
*
* Revision 1.1.1.1  2002/07/24 15:56:25  rdm
* initial import into CVS
*
* Revision 1.1.1.1  2002/06/16 15:18:39  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:17  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:56  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 02/07/94  18.24.47  by  S.Giani
*-- Author :
      SUBROUTINE G3SROTM(NMAT,THETA1,PHI1,THETA2,PHI2,THETA3,PHI3)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       STORE ROTATION MATRICES                                  *
C.    *                                                                *
C.    *    ==>Called by : <USER>                                       *
C.    *         Author  R.Brun  *********                              *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gcunit.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcnum.inc"
      COMMON / FIXIT / JR
      DIMENSION ANGLES(6),IP(3)
      SAVE SINMIN
#if defined(CERNLIB_SINGLE)
      DIMENSION ROTMAT(9)
      PARAMETER (ONE=1.0,ZERO=0.0)
      DATA SINMIN/1.00E-5/
#endif
#if !defined(CERNLIB_SINGLE)
      DOUBLE PRECISION ROTMAT(9),ONE,ZERO
      DOUBLE PRECISION PROD1,PROD2,PROD3,HMOD,SINMIN
      PARAMETER (ONE=1.D0,ZERO=0.D0)
      DATA SINMIN/1.00D-5/
#endif
C.
C.    ------------------------------------------------------------------
C.
      IF(NMAT.LE.0)GO TO 999
      IF(JROTM.LE.0)CALL MZBOOK(IXCONS,JROTM,JROTM,1,'ROTM',NROTM,NROTM,
     +0,3,0)
      IF(NMAT.GT.IQ(JROTM-2)) THEN
         NPUSH=NMAT-IQ(JROTM-2)+50
         CALL MZPUSH(IXCONS,JROTM,NPUSH,0,'I')
         NROTM=IQ(JROTM-2)
         JR1=0
      ELSE
         JR1=LQ(JROTM-NMAT)
         IF(JR1.GT.0)THEN
            WRITE(CHMAIL,1000)
            CALL GMAIL(1,0)
            CALL G3PROTM(NMAT)
            CALL MZDROP(IXCONS,LQ(JROTM-NMAT),' ')
         ENDIF
      ENDIF
      CALL MZBOOK(IXCONS,JR,JROTM,-NMAT,'ROTM',0,0,16,3,0)
C
      Q(JR + 11) = THETA1
      Q(JR + 12) = PHI1
      Q(JR + 13) = THETA2
      Q(JR + 14) = PHI2
      Q(JR + 15) = THETA3
      Q(JR + 16) = PHI3
C
      DO  10 N = 1,3
        THERAD = Q(JR+ 9+2*N)*DEGRAD
        PHIRAD = Q(JR+10+2*N)*DEGRAD
        SINTHE = SIN(THERAD)
        Q(JR+3*N-2) = SINTHE * COS(PHIRAD)
        Q(JR+3*N-1) = SINTHE * SIN(PHIRAD)
        Q(JR+3*N  ) = COS(THERAD)
        CALL VUNIT (Q(JR+3*N-2),Q(JR+3*N-2),3)
  10  CONTINUE
C.
C.---       Test orthonormality
      DO 20 J=1,9
        ROTMAT(J)=Q(JR+J)
  20  CONTINUE
      PROD2=ZERO
C.
C.               X - Y
      PROD1=
     +ROTMAT(1)*ROTMAT(4)+ROTMAT(2)*ROTMAT(5)+ROTMAT(3)*ROTMAT(6)
      IF(ABS(PROD1).GT.SINMIN) GO TO 30
C.
C.               X - Z
      PROD2=
     +ROTMAT(1)*ROTMAT(7)+ROTMAT(2)*ROTMAT(8)+ROTMAT(3)*ROTMAT(9)
      IF(ABS(PROD2).GT.SINMIN) GO TO 30
C.
C.               Y - Z
      PROD3=
     +ROTMAT(7)*ROTMAT(4)+ROTMAT(8)*ROTMAT(5)+ROTMAT(9)*ROTMAT(6)
      IF(ABS(PROD3).LE.SINMIN) GO TO 110
  30  CONTINUE
C.
C.---       Orthonormalization needed
C.
C.          Assume X correct
      HMOD=ZERO
      DO 40 J=4,6
        ROTMAT(J)=ROTMAT(J)-ROTMAT(J-3)*PROD1
        HMOD=HMOD+ROTMAT(J)*ROTMAT(J)
  40  CONTINUE
      HMOD=ONE/SQRT(HMOD)
      DO 50 J=4,6
        ROTMAT(J)=ROTMAT(J)*HMOD
  50  CONTINUE
C.
C.          Y done, do Z
C.
*      IF(PROD2.EQ.ZERO) THEN
*        PROD2=
*     +  ROTMAT(1)*ROTMAT(7)+ROTMAT(2)*ROTMAT(8)+ROTMAT(3)*ROTMAT(9)
*      ENDIF
*      PROD3=
*     +ROTMAT(4)*ROTMAT(7)+ROTMAT(5)*ROTMAT(8)+ROTMAT(6)*ROTMAT(9)
*      HMOD = ZERO
*      DO 60 J=1,3
*        ROTMAT(J+6)=ROTMAT(J+6)-ROTMAT(J+3)*PROD3-ROTMAT(J)*PROD2
*        HMOD = HMOD+ROTMAT(J)*ROTMAT(J)
*  60  CONTINUE
* == AV ==
      ROTMAT(7) = ROTMAT(2)*ROTMAT(6) - ROTMAT(3)*ROTMAT(5)
      ROTMAT(8) = ROTMAT(3)*ROTMAT(4) - ROTMAT(1)*ROTMAT(6)
      ROTMAT(9) = ROTMAT(1)*ROTMAT(5) - ROTMAT(2)*ROTMAT(4)
      HMOD = ZERO + ROTMAT(7)*ROTMAT(7) + ROTMAT(8)*ROTMAT(8)
     &            + ROTMAT(9)*ROTMAT(9)
* == AV ==
      HMOD = ONE/SQRT(HMOD)
      DO 70 J=7,9
        ROTMAT(J) = ROTMAT(J)*HMOD
  70  CONTINUE
C.
C.          Put back the matrix in place
      DO 80 J=1,9
        Q(JR+J)=ROTMAT(J)
  80  CONTINUE
C.
C.          Now recompute the angles
      DO 90 J=1,3
        ANGLES(J*2-1) = ACOS(MAX(-ONE,MIN(ONE,ROTMAT(J*3))))*RADDEG
        ANGLES(J*2)   = ZERO
        IF(ROTMAT(J*3-1).NE.ZERO) THEN
          ANGLES(J*2) = ATAN2(ROTMAT(J*3-1),ROTMAT(J*3-2))*RADDEG
          IF(ANGLES(2*J).LT.0.0) ANGLES(2*J) = ANGLES(2*J)+360.0
        ENDIF
  90  CONTINUE
      WRITE(CHMAIL,2000) NMAT
      CALL GMAIL(1,2)
      WRITE(CHMAIL,2001) (Q(JR+10+J),J=1,6)
      CALL GMAIL(0,0)
      WRITE(CHMAIL,2002) ANGLES
      CALL GMAIL(0,1)
C.
C.          Put back the angles in place
      DO 100 J=1,6
        Q(JR+10+J) = ANGLES(J)
 100  CONTINUE
C.
C.---       Orthonormalization ended
 110  CONTINUE
C
      DO 130 J = 1,3
        IP(J)  = 3
        JJR=JR+J*3-3
C
        DO 120 I = 1,3
          IF(ABS(Q(JJR+I)).LT.0.99999) GO TO 120
C
          IP(J)  = I + 3
          IF(Q(JJR+I).GE.0.) GO TO 130
C
          IP(J)  = 3 - I
          GO TO 130
C
 120    CONTINUE
 130  CONTINUE
C
      Q(JR + 10) = IP(1) + 10* IP(2) + 100* IP(3)
C
      IF(JR1.GT.0) THEN
         CALL G3PROTM(-NMAT)
      ENDIF
C
1000  FORMAT(' *** GSROTM ***: Warning, rotation matrix redefinition:')
2000  FORMAT(' *** GSROTM ***: ',
     +       'Parameters of matrix no. ',I4,' changed:')
2001  FORMAT(' Old values: ',6(F14.7,3X))
2002  FORMAT(' New values: ',6(F14.7,3X))
 999  RETURN
      END
